(*
#use "MLlab.ml";; (* inclus le code de MLlab.ml *)

#use "topfind";; (* utilisation de findlib *)
#require "lacaml";; (* charge le package 'lacaml' *)
#load "MLlab.cmo";;
#load "Lacaml_ext.cmo";;
#load "Lacaml_L1.cmo";;
#load "projection.cmo";;
*)
open Lacaml.D;;

(***** test de la projection L1 ******)
(* on doit trouver [3; 4; -13; 0; 0] *)
let n = 5;;
let c = Lacaml.D.Vec.make0 5;;
begin
c.{1} <- 28.0;
c.{2} <- 29.0;
c.{3} <- -38.0;
c.{4} <- 21.0;
c.{5} <- -4.0;
end;;
let tau = 20.0;;
let x = Projection.projector tau c;;

(***** autre test de la projection L1 *****)
let c =  Lacaml_ext.D.Vec.randn 10;;
let tau = 2.0;;
let x = Projection.projector tau c;;
let nrm1_x = Lacaml_ext.D.Vec.nrm1 x;;


(* build sparse spikes vector
   example :
   let x0 = make_sparse_spikes 10 4 (fun () -> MLlab.Random.randn 0.0 1.0);;
*)
let make_sparse_spikes ?rnd_state n k amplitude_law =
  let vec = Vec.make0 n in
  let state = 
    match rnd_state with
        | None -> Random.get_state ()
        | Some state -> state in
  let p = MLlab.Random.randperm ~rnd_state:state n in
  for j = 1 to k
  do
    vec.{p.(j-1)} <- amplitude_law ();
  done;
  if rnd_state = None then Random.set_state state;
  vec;;

(* build random encoding matrix *)
let make_random_encoding m n =
  let a = Lacaml_ext.D.Mat.randn n m in (* noter la transposition *)
  let tau = geqrf a in
  begin
    orgqr ~tau:tau ~k:m a;
    Mat.transpose a;
  end;;

(* m rows, n cols, k nonzeros *)
let m = 4 and n = 6 and k = 2 in
let a = make_random_encoding m n in
let x0 = make_sparse_spikes n k (fun () -> MLlab.Random.randn 0.0 1.0) in
Format.printf "@[<2>A =@\n@\n@[%a@]@]@\n@\n" pp_mat a;;
